### Replication Code for:
# "Preferences over Foreign Migration: Testing Existing Explanations in the Gulf"
# World Politics 2022
# All figures and tables generated in this file use data from:
# World Values Survey Round 7 and
# UN population division global migrant stock (2019)

### Author Erin York
### File updated February 3, 2022
### File created using R version 4.1.1 (2021-08-10) -- "Kick Things"


# 1 Load data and required packages ---------------------------------------

rm(list = ls())

library(tidyverse)
library(sandwich)

# set working directory as appropriate
setwd("YOUR_FILEPATH_HERE")

# load in WVS wave 7
load("WVS_Cross-National_Wave_7_R_v2_0.rdata")

wvs<- `WVS_Cross-National_Wave_7_v2_0`

# load in UN migrant stock data
undat<- read.csv("un_migrantpop.csv")


# 2 Merge and recode datasets ---------------------------------------------


merge_wvs<- wvs %>% 
  rename(Code = B_COUNTRY) %>%
  # add NAs where relevant in WVS data
  mutate(
    # migrant unemployment attitudes  
      migrants_increaseunemploy = ifelse(Q128 < 0, NA, Q128),
         migrants_against = ifelse(Q130 < 0, NA, Q130),
         # demographics
         female = ifelse(Q260 < 0, NA, ifelse(Q260 == 2, 1, 0)),
         age = ifelse(Q262 < 0, NA, Q262),
         immigrant = ifelse(Q263 < 0, NA, ifelse(Q263 == 2, 1, 0)),
         education = (ifelse(Q275 < 0, NA, Q275)),
         employ = ifelse(Q279 <= 0, NA, ifelse(Q279 >3, 0, 1)),
         unemploy = ifelse(Q279 ==7, 1, 0)) %>%
  left_join(undat) %>%
  # remove immigrants from sample
  filter(immigrant == 0)


# 3 Generate figures and analysis -----------------------------------------

####
## Figure A3 ##


undat2<- undat %>% 
  select(Code, migrant_2019) %>% unique() %>%
  # add facet label to UN data
  mutate(type = "All Countries")

# Figure output
merge_wvs %>% 
  select(Code, migrant_2019) %>% 
  unique() %>%
  # add facet label to WVS data
  mutate(type = "WVS Sample") %>%
  # join UN and WVS data
  bind_rows(undat2) %>%
  # generate histogram
  ggplot(aes(x = migrant_2019)) +
  geom_histogram() + 
  facet_grid(type~., scales = "free_y") + 
  theme_bw() +
  ylab("") + 
  xlab("Foreign % of Population (2019)")

####

####
## Table A10 ##

# model 1
against1<- (lm(migrants_against ~ female + age  + education + (migrant_2019) +
          factor(regionWB) + log(popWB2019), merge_wvs))

# model 2
against2<- (lm(migrants_against ~ female + age  + education + log(migrant_2019) +
           factor(regionWB) + log(popWB2019), merge_wvs))

# model 3
unemploy1<- (lm(migrants_increaseunemploy ~ female + age + education + (migrant_2019) +unemploy+
          factor(regionWB) + log(popWB2019), merge_wvs))

# model 4
unemploy2<- (lm(migrants_increaseunemploy ~ female + age  + education + log(migrant_2019) +unemploy+
           factor(regionWB) + log(popWB2019), merge_wvs))

# model 5
unemploy_int1<- (lm(migrants_increaseunemploy ~ female + age + education + (migrant_2019)*unemploy +
          factor(regionWB) + log(popWB2019), merge_wvs))

# model 6
unemploy_int2<- (lm(migrants_increaseunemploy ~ female + age  + education + log(migrant_2019)*unemploy +
           factor(regionWB) + log(popWB2019), merge_wvs))


# generate Table A10
stargazer::stargazer(against1, against2, unemploy1, unemploy2, unemploy_int1, unemploy_int2,
                     # SEs clustered at country level (n clust = 51)
                     se = list(sqrt(diag(vcovCL(against1, merge_wvs$Code))),
                               sqrt(diag(vcovCL(against2, merge_wvs$Code))),
                               sqrt(diag(vcovCL(unemploy1, merge_wvs$Code))),
                               sqrt(diag(vcovCL(unemploy2, merge_wvs$Code))),
                               sqrt(diag(vcovCL(unemploy_int1, merge_wvs$Code))),
                               sqrt(diag(vcovCL(unemploy_int2, merge_wvs$Code)))),
                     keep = c("migrant_", "unem*"),
                     covariate.labels = c("Migrant \\% of Pop", "Log Migrant \\%",
                                          "Unemployed",
                                          "Migrant \\% x Unemployed", "Log Migrant \\% x Unemployed"),
                     omit.stat = c("rsq", "f", "ser"),
                     dep.var.labels  =c("Restrictions", "Increase Unemployment", "Increase Conflict"),
                     
                     add.lines = list(c("Controls", rep("\\checkmark", 6))),
                     float = F)

####


